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We discuss numerical solution of Altarelli-Parisi equations in a Laguerre-polynomial 
method and in a brute-force method. In the Laguerre method, we get good accuracy 
by taking about twenty Laguerre polynomials in the flavor-nonsinglet case. Excellent 
evolution results are obtained in the singlet case by taking only ten Laguerre polyno- 
mials. The accuracy becomes slightly worse in the small and large x regions, especially 
in the nonsinglet case. These problems could be implemented by using the brute-force 
method; however, running CPU time could be significantly longer than the one in the 
Laguerre method. 

1. Introduction 

In order to investigate nucleon substructure, we scatter high-energy electrons or 
muons from the proton target. Its cross section is expressed in term of two unpolarized 
structure functions and F2, which depend on two kinematical variable x and Q^. 
If the nucleon consists of free partons, the structure functions are independent of the 
variable . This assumption is called Bjorken scaling hypothesis. However, it is well 
known that the structure functions do depend on Q"^ experimentally and theoretically. 
This fact is called "scaling violation" , which is important in testing perturbative QCD. 
The scaling violation results from the fact that high-energy photons probe minute qq 
clouds around a quark. 

An intuitive way to describe the phenomena is to use the Altarelli-Parisi equation. 
The fiavor-nonsinglet Altarelli-Parisi equation is given by 



where q^si^^Q^) ^ nonsinglet quark distribution, asiQ"^) is the running coupling 
constant, and PNsi^) is a nonsinglet splitting function. The above equation describes 
the process that a quark with the nucleon's momentum fraction y radiates a gluon and 
it becomes a quark with the momentum fraction x. The splitting function Pns{z) de- 
termines the probability for a quark radiating a gluon such that the quark momentum 
is reduced by the fraction z. Next-to-leading-order corrections to the Altarelli-Parisi 
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equation are included in the coupling constant as and in the splitting function Pns- 
The leading order and the next-to-leading order are abbreviated to LO and NLO re- 
spectively throughout this paper. 

The singlet Altarelli-Parisi equations are coupled integrodifferential equations: 
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where ^^(x, Q"^) is the singlet quark distribution = ^(g^ + Qi) and G{x, Q^) is the 
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gluon distribution. Pgg{x), Pgcix), Pcqix), and Pgg{x) are splitting functions. 

Because the evolution equations are very important for testing perturbative 
QCD and are often used theoretically and experimentally, it is worth while having a 
computer program of solving it accurately without consuming much computing time. 
Our research purpose is to create a useful computer program for numerical solution 
of the Altarelli-Parisi equations. There are several methods in solving the equations; 
however, a Laguerre-polynomial method is considered to be very effective among them 
0. In Ref. 0, we provide the Laguerre method program LAG2NS, which deals with 
the Q2 evolution of a nonsinglet distribution in the next-to-leading order. In this 
paper, we show Laguerre-method results in Ref. and compare them with results in 
a brute-force method @]. 

In section 2, we explain the Laguerre method for solving the evolution equations 
with NLO corrections. The brute-force method is discussed in section 3 and conclusions 
are given in section 4. 



2. Laguerre-polynomial method 

2.1 Nonsinglet distribution 

The Laguerre-polynomial method is considered to be effective in convergence and 
in computing time. This is because the first several polynomials resemble parton distri- 
butions. Numerical solution for the evolution equation including the NLO corrections 
is discussed in the following. We refer the reader to the papers in Refs. Q for more 
complete account of the Laguerre method. The NLO running coupling constant is used 
in this paper unless we specify LO. The renormalization scheme is MS. 

We first define the variable t and the evolution operator E(x,t) hj t = ——In 

and 

q^,{x,t)= ^ E{-,t)q^,{y,t = 0) , (2.1) 
Jx V V 
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where q is given by q{x) = xq{x). The function E{x, t) describes the evolution from 
the initial distribution at Qg (t = 0) to the one at Q^. Substituting Eq. (2.1) into the 
Altarelli-Parisi equation in Eq. (1.1), we find that the evolution operator satisfies 



— E{x,t) = / — 
dt Jx y 



y 



2tt 



y 



E{y,t) 



(2.2) 



where P{x) = xP{x) and R{x) are given by R{x) = Pj^li^) ^^Nsi^)- ^Nsi^) ^^"^ 

P^l{x) are LO and NLO splitting functions. The function E{x,t) is split into LO 
and NLO contributions by E{x,t) = E'-°\x,t) + ^^^£^(^)(x, t). Substituting this into 

Eq. (2.2) and considering that the LO evolution operator E^^^{x,t) satisfies the LO 
Altarelli-Parisi equation, we find that the LO and NLO evolution operators are not 
independent and they are related by 
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The variable x' is introduced by x' = — In x (0 < x < 1) in order to use the Laguerre 
polynomials, which are defined in the region < x' < oo. Using the variable x', we ex- 
pand the functions in terms of the Laguerre polynomials: f{x = ) = fnLnjx'). 

n 

Then, the problem is to obtain the expansion coefficients En in terms of the coeffi- 
cients Pn of the splitting function. The LO evolution operator is given in Ref. [ffl] 
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by E^^\t) = e-^o* J2 -j^_^n, ^^^^e is calculated by the recurrence relation B, 
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i=k j=l 

The NLO evolution operator is given by using the calculated E^^^t) and Eq. (2.3): 
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where Ri are the Laguerre coefficients of R{x), and is defined as 
these results, we finally obtain 
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0. Using 

(2.5) 



«s(0) 



with En{t) = E^^\t) + ^£;(i)(t). In this way, the original integrodifferential equa- 

tion is reduced to a sum of finite number of Laguerre expansion coefficients and the 
Laguerre polynomials. 
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We have discussed the evolution of a nonsinglet quark distribution. evolution 
of a nonsinglet structure function is calculated in a similar way. The only modification 
is to take into account NLO corrections from the coefficient function. A nonsinglet 
structure function F^si^, Q"^) is expressed as a convolution of the corresponding non- 
singlet quark distribution q^^g^x^Q'^) and the coefficient function CNs{x,as) as 



y y 



(2.6) 



Cns is given by C^six^as) = 5{1 - x) + -^B--[x 
tion. We show examples of our evolution results 

The initial distribution is given by 
the HMRS-B valence-quark distribution 
at Q2=4 GeV^ [§. This distribution 
is evolved to the one at Q^=200 GeV^. 
Evolved distributions are shown in Fig. 

1. We find that accurate results are ob- 
tained by taking about twenty Laguerre 
polynomials; however, the accuracy be- 
comes slightly worse in the very small x 
region {x < 0.001) and in the very large x 
region {x > 0.9). 

Our evolution results are compared 
with the CDHSW neutrino data in Fig. 

2. The CCFR ^3(0;) distribution at Q^=3 
GeV^ ^ is chosen as the initial distribu- 
tion in both LO and NLO cases. We em- 
ploy A=0.21 GeV and four flavors for cal- 
culating the evolution. The NLO contri- 
butions are important at small (~ 1 
GeV^) and at small shown in Fig. 
2. Furthermore, the flgure indicates that 
our variations are consistent with the 
CDHSW experimental data. 



, where B (x) is the NLO correc- 
in the following. 



EPS File figl.eps.ps not 
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Fig. 1 evolution of x{uy + d^)- 
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Fig. 2 Q'^ evolution of F3 
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2.2 Singlet distribution 

The nonsinglet evolution is rather simple because gluons do not participate in the 
evolution process. However, the singlet part is more complicated because the Altarelli- 
Parisi equations are coupled integrodifferential equations in Eq. (1.2). We do not 
discuss formalism of singlet solution, so interested reader may read the papers in Ref. 
[Q. Even though the analysis is slightly complicated, the essential procedure is the same 
as the one in the previous subsection. Our investigation for the singlet part is still in 
progress 0], so the results in this subsection should be considered as preliminary ones. 
Calculating Laguerre-expansion coefficients for the splitting functions in Eq. (1.2) and 
for initial distributions, we obtain a solution similar to Eq. (2.5). 

The initial singlet distribution qs{x) and gluon distribution G{x) are chosen as the 
HMRS-B distributions at Q^=A GeV^. They are evolved to the ones at ^^=200 GeV^. 
Three evolution results are shown by taking five, ten, and twenty Laguerre polynomials 
in Figs. 3 and 4. It is rather surprising to find excellent accuracy even at very small x 
as small as 10~^. The accuracy is much better than the nonsinglet one by taking only 
ten Laguerre polynomials. 



EPS File fig3.eps. ps not 
found 



Fig. 3 evolution of xqs{x) 



EPS File fig4.eps. ps not 
found 



Fig. 4 Q'^ evolution of xG{x). 



Typical running time is a few seconds on SUN-IPX or on VAX-4000/500 in obtaining 
a x-distribution curve in Fig. 1, 3, or 4. Considering the excellent accuracy and 
the short running time, we conclude that the Laguerre method is very effective for 
numerical solution of the Altarelli-Parisi equations. 
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3. Brute-force method 



The Laguerre method is practically useful; however, it does not produce excellent 
results at small x and at large x in the nonsinglet case. With development of high- 
energy accelerators, the small x region becomes more and more interesting. For those 
who are serious about the small x and large x regions, we discuss an alternative method, 
the brute- force method in Ref. 0], for obtaining accurate numerical solution. 

The brute-force method is perhaps the simplest one in solving the integrodifferential 
equation. We divide the variables t and x into small steps, then the derivative by t is 
defined as 

d q{xi,tj+i) - q{x^,tj) 

^ — , [6.L) 

and the integral over y is by 



' ^ Pi-) q{y,t) ^ f: ^ P{^) q{xk,t,) . (3.2) 

y y k^i) Xk xk 



If the number A^^, is increased, we expect to obtain accurate numerical results. However, 
numerical convergence is not good in the small x region even if we take large A^^^. So we 
decide to divide x into equal steps in the logarithmic scale. Then, our evolution results 
converge by taking large A^^; ~ 200. Evolution results of the HMRS-B valence-quark 
distribution are shown in Fig. 5 by taking A^j;=50, 200, 500 with fixed number of points 
A'^t=100 for the variable t. The accuracy is excellent in the nonsinglet distribution, so 
the brute-force method could be used in the small x region instead of the Laguerre 
method. However, it takes significant amount of running CPU time as shown in Fig. 
6. 



EPS File fig5.eps. ps not 
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Fig. 5 evolution of x{uv + d^)- 



Fig. 6 Computing time. 
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4. Conclusions 



We investigated numerical solution of the Altarelli-Parisi equations with next-to- 
leading-order corrections in the Laguerre-polynomial and the brute-force methods. Ex- 
cellent convergence is obtained in the singlet case by using the Laguerre method with 
merely ten Laguerre polynomials. The accuracy is slightly worse in the nonsinglet case 
especially at small and large x. However, accurate results are obtained in these regions 
by using the brute-force method. 
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